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We also describe an implementation of the algorithm in a Soften Lagrangian Hy- 
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well as a realistic simulation of cosmological reionization. 
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1 Introduction 



Numerous questions in pliysical cosmology and galaxy formation require a 
detailed understanding of the physics of radiative transfer. Of particular in- 
terest are problems related to the reionization of the intergalactic medium 
(Gnedin & Ostriker, 1997; Haiman et al., 2000; Barkana & Loeb, 2001, and 
references herein), absorption line signatures of high redshift structures (Abel 
& Mo, 1998; Kepner et al., 1999), and hydrodynamic and thermal effects of 
UV emitting central sources on galaxy formation (e.g. Haehnelt, 1995; Silk & 
Rees, 1998). 

Several attempts have been made to develop numerical algorithms for solving 
the radiative transfer equation in astrophysics in general, and in cosmology 
in particular. Abel et al. (1999) described a ray tracing approach for model- 
ing radiative transfer around point sources. A related ray tracing approach 
was introduced by Razoumov & Scott (1999) which can also handle diffuse 
radiation fields but finite computer memory limits the angular resolution and 
therefore only poorly captures ionization fronts around point sources on rea- 
sonably large grids. Ciardi et al. (2001) discussed a Monte Carlo approach 
to sample angles, energies and the time evolution of the radiative transfer 
problem. A ray tracing scheme for use with Smooth Particle Hydrodynam- 
ics was introduced by Kessel-Deynet & Burkert (2000). All these approaches 
however suffer from the excessive operation count and are costly to implement 
in reahstic simulations Sokasian et al. (e.g. 2001). 

In order to implement a radiative transfer scheme into modern numerical 
simulations, one needs to be able to solve the radiative transfer equation in 
0{N) operations (where is the number of baryonic resolution elements, and 
we ignore a possible logarithmic multiplier, so that 0[N log N] still counts as 
0[A^]). However, the radiative transfer equation is six-dimensional in general 
and its direct solution requires of the order of 0{N^^^ x N^) operations, where 
is the number of frequency bins (typically must be about 200-300). Various 
formal solutions of the transfer equation (like an attenuation equation) can be 
reduced to a four dimensional (spatial plus frequency dimensions) form, but 
then the calculation of the pairwise optical depth requires an operation count 
of the order of 0{N x Ns x Ny)^ where Ns is the number of sources, which for 
an arbitrary source function can be equal to A^. For some restricted case, like 
reionization by quasars, when the number of sources Ns is small enough, the 
whole solution can indeed be computed directly in 0{N) operations (Abel et 
al., 1999). However, for more realistic cases the operation count increases to 
0{N'^). In addition, in cosmological applications calculation of the background 
radiation field would still require N^^^ operations. 

Thus, one has to use approximations to try to decrease the operation count. 
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The "Local Optical Depth" (LOD) approximation for the cosmological radia- 
tive transfer was first introduced in Gnedin & Ostriker (1997) for the homo- 
geneous distribution of sources, and later was expanded on the case with an 
arbitrary source function in Gnedin (2000). This approximation is however 
highly approximate and uncontrolled. 

Another possible approach is to use a so-called Variable Eddington Tensor 
approximation by solving moment equations for the radiative field. Such an 
approach allows to ensure exact conservation of photon numbers and flux, 
and thus would automatically give correct speeds for the ionization fronts. In 
this paper we describe such an approach whose main feature is the Variable 
Eddington Tensor approximation with the tensor calculated in the optically 
thin regime. 

Norman et al. (1998) first suggested to solve the moment equations for the 
simulation of three dimensional radiative transfer in cosmological hydrody- 
namics. They suggested to split the transfer of radiation of point sources from 
the transfer of diffuse emitters. The method presented here shares some of 
their ideas, but works with different variables. We reformulate the equations 
of radiative transfer by decomposing the radiation field into a uniform compo- 
nent and a perturbation term. It is this new mathematical framework that we 
will present in the following section. In Section 3 we discuss our treatment of 
multiple frequency bins and also introduce a new "Reduced Speed of Light" 
approximation which allows us to overcome numerous numerical problems. We 
then go on in Section 4 to introduce particular details of our implementation 
of the algorithm. Section 5 is devoted to various one and three dimensional 
test cases that exemplify the strengths and weaknesses of our approach. We 
also present in that section a realistic cosmological simulation based on our 
new scheme, and compare our scheme with the LOD approach. We conclude 
in section 6 and give a discussion of various performance issues. 



2 Formalism 



2.1 Basic equations 



The evolution of the specific intensity J^, in the expanding universe is given 
by the following equation: 



dt 



B f B J \ 



3 



Here ,x* are the comoving coordinates, H is the Hubble constant, k^, is the 
absorption coefficient, is the source function, and = cn^ /a, where is 
the unit vector in the direction of photon propagation. 

Equation (1) inchides all relativistic effects and therefore is quite complex. But 
in the same spirit as equations for particle motions in the expanding universe 
are split into a uniform relativistic background and Newtonian perturbations, 
we can split equation (1) into a background and perturbations (Gnedin & 
Ostriker, 1997). 

First, we define the mean specific intensity, 

J,{t) = {J,{t,x,n)), (2) 



where the averaging operator () acting on a function f{x,n) of position and 
direction is defined as: 

{fix, n)) =lim^^l £x I dnfix, n). (3) 

V 



The mean specific intensity ( "background" ) satisfied the following equation: 
^-H l^u^ - 3X) = -hi + S^, (4) 

where, hy definition, 



and 



r {ki/Jf) 

= = . 



It is important to emphasize here that we use the overbar symbol merely as 
a part of notation and not to denote the space average; in particular, h is 
not a space average of h since it is weighted by the local value of the specific 
intensity J^. 

We can now introduce a new quantity, which we call relative (to the mean) 
specific intensity /^(t, x, n), which quantifies the difference between the local 
value of the specific intensity and the mean. 



Jv — fvJvi (5) 



4 



so that (fi,) = 1. In the dynamic analogy, fi, corresponds to 1 + 5, where 5 is 
the cosmic overdensity. 



It is now straightforward to derive the following equation for f„: 



This equation is of course not simpler than the original equation (1). But, if 
we restrict ourselves to scales significantly smaller than the horizon size, and 
matter velocities much smaller than the speed of light (Newtonian limit), then 
we can simplify equation (6) substantially by noting that the relative specific 
intensity does not change substantially and the universe does not expand 
substantially over the period of time a photon needs to travel over the scale 
under consideration (for example, a size of a computational box). Therefore, 
we can ignore the derivative with respect to the frequency. Strictly speaking, 
we should also discard the time derivative, as it has the same power of v/c as 
the frequency derivative. However, we will retain it for a moment, keeping in 
mind that it is small - as we explain below, the time derivative turns out to 
be quite useful for the numerical implementation of our scheme. 

Under these assumptions, equation (6) reduces to the following simple equa- 
tion: 

adf^ w ^ / 



where 



1^1/ — ki/} 
c 



Kn — Kn Ku H =- (8) 



and 



a 

c J„' 



Equations (4) and (7) thus separate the average cosmological evolution of 
specific intensity (including redshift and effects of "distant" sources - sources 
outside the simulation box) from the local evolution due to inhomogeneous 
absorption and sources inside the box. 
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2.2 Moments of the transfer equation 



The transfer equation (7) is a six-dimensional equation in the phase space, 
and as such offers a considerable difficulty to solve numerically. Pursuing the 
analogy with dynamics further, we can define moments of the radiation field 
similarly to how density, momentum, pressure etc of a medium are defined as 
moments of the distribution function. Specifically, we define the (relative to 
the mean) radiation energy density E^, 

E^{t,x) = j dD.f^{t,x,n), 

the flux F^, 

F^{t,x) = ^JdQn'f,{t,x,n), 

and the so-called Eddington tensor h]^ , 

E^{t,x)h}^{t,x) = j dfln'nf^{t,x,n). 

In particular, 

Tr/it^ = l and (E^) ^ 1. 

Taking moments of equation (7), we can derive equations for the moments of 
the radiation field (assuming that the source function is isotropic): 

adE^ dFi . ^ , 

and 

adF^ d ■ 

It is easy to see that equations (9) and (10) are in a conservative form, and so 
represent the conservation of the number of photons and the flux respectively. 

Equations (9) and (10) constitute a system of two partial differential equations. 
The system is however not closed since the Eddington tensor is not specified by 
these two equations. However, if the Eddington tensor is given as a function 
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of position and time, the system becomes closed. In the following section 
we assume that the Eddington tensor is known and describe the numerical 
implementation of an algorithm based on equations (9) and (10). We then 
discuss how to define the Eddington tensor, thus completing the numerical 
scheme. 



3 Specific implementation 



3.1 The "Reduced Speed of Light" approximation 



In this section we omit the frequency dependence for clarity. Equations (9) 
and (10) can be combined into one equation preserving the accuracy of the 
Newtonian approximation. Namely, from (10) we have: 



k \ c dt dx 



Substituting this into (9), we obtain: 

adE d (ldEh'i\ a d fldF^\ 

__ + ,^_^__ ___ . (11) 



The right hand side can be expanded into two terms, 

ad (ldF^\_ akjdF^ a d dF^ 
c dx^ \k dt J CK^ dt ck dt dx^ 



Both of these terms are of the order of 1/c, and thus can be dropped in the 
Newtonian limit. We therefore neglect the first term, but retain the second 
one, keeping in mind that it is of the higher order in 1/c than the leading 
terms in equation (11). That leaves the right hand side of equation (11) in the 
following form: 

a^d_ f_ad^_^^^\ 
ckdt \ c dt / ' 



of which we retain only one term. 



adE 
'c~dt' 



7 



with the same 1/c accuracy. This reduces the original equation (11) to the 
following final form: 



This one equation is equivalent to the original two equations (9) and (feq) 
in the Newtonian limit, and the term with the time derivative is small and 
formally of the order of 1/c. 

The physical reason for retaining the time derivative and throwing away other 
terms which are formally of the same order is that at the ionization front E 
can change very rapidly making this term significantly larger than the terms 
we neglected, even if formally they have the same 1/c order. 

In the Newtonian limit the speed of light is infinite, which in practice means 
that we ignore all the terms of the order of v/c. Thus, it actually does not 
matter what the specific value of the speed of light is as long as f/c <^ 1, and 
we can replace the true speed of light c in equation (12) with some other quan- 
tity c without violating the validity of the Newtonian limit as long as v/c <^ 1. 
For example, for modeling cosmological reionization we can adopt a value for 
c as low as 1,000 km/ s, because typical gas velocities during reionization do 
not exceed 100 km/ s. 

This trick of "reducing the speed of hght" (which we call the Reduced Speed 
of Light, or RSL, approximation) allows us to achieve a very important goal. 
Namely, if we ignore the time derivative term in equation (12), the equation 
turns elliptic, which is difficult to solve for an arbitrary k. Thus, keeping the 
time derivative in equation (12) allows us to solve this equation as a dynamical 
equation and use simple stable explicit schemes. 

There is however one very serious limitation - the time-step in any explicit 
scheme will be limited by an appropriate Courant condition to a value which is 
about cs/c times smaller than the hydrodynamic time-step (here cs is the gas 
sound speed). For the intergalactic gas at 10, 000 K this ratio can be as small as 
1/30, 000. Clearly, this will make such a scheme impractical, as it is impossible 
to make 30,000 radiative transfer time steps for every hydrodynamic time step 
in a realistic simulation. 

The RSL approximation allows to reduce this ratio dramatically to less than 
100, and in our tests described below we even find that taking c ~ does 
not lead to large errors. 




(12) 
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Thus, in the RSL approximation equation (11) reduces to the following form: 



where c < c. 

It is important to note here that the reduction in the speed of light does not 
limit the propagation speed of an ionization front to c. Equation (13) is a 
diffusion equation, which has an infinite speed of signal propagation. In other 
words, the RSL approximation puts no limit on the speed of ionization fronts 
and thus can sometimes lead to unphysical, faster than the speed of light 
propagation of the ionization fronts (Abel et al., 1999). In practice, we put a 
limit on the maximum possible speed of an ionization front by limiting the 
number of time-steps. This is discussed in greater detail in §3.4. 

3.2 Multi-frequency transfer 

Equation (13) is a diffusion equation, and as such is not more difficult to solve 
than equations of hydrodynamics (in fact, it is much easier to solve). However, 
we have omitted the frequency dependence in the derivations of the previous 
section. Tests indicate that in order to calculate the various ionization and 
heating rates with sufficient precision on a logarithmically spaced mesh in 
frequency space, at least 20 points per e-folding are required, which translates 
into at least 200 frequency bins. Solving equation (13) 200 times at each 
hydrodynamic time step would be too prohibitive with the currently available 
computer resources. 

We therefore need to introduce further approximation in order to break this 
unfavorable scaling. In order to do that, it is instructive to consider the op- 
tically thin case first. In this regime the angle averaged solution to equation 
(1) is (omitting time dependence for clarity): 



Both, J,^ and S^, may have very complicated frequency dependence, which 
would require a large number of frequency bins. However, if we assume that 
all sources have the same frequency dependence, i.e. 



where depends only on frequency, and depends only on time and space 




(13) 




S„{t,x) = g„p^{t,x), 
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(say, the mass density of stars), then in the optically thin regime J,^ as a 
function of frequency separates into two terms: 



(14) 



where Ei — 1 and 



D2 



a 



I 




(15) 



and both Ei and D2 are independent of frequency. Thus, in this regime the 
frequency dependence is factorized in two one dimensional functions Ji, and 
Qi,, and all relevant ionization and heating rates need to be computed only 
once and then simply rescaled by frequency-independent factors Ei and D2. 

This procedure is equivalent to separating all sources into "distant" ones, 
which include full redshift dependence and cosmological terms, as well as dif- 
fuse radiation, and whose frequency dependence is described by J^,, and into 
"near" sources which are just Newtonian perturbations on the GR background 
with the frequency dependence given by Qi,. 

In a general case, let's try to use ansatz (14). In that case Ei and D2 are 
no longer frequency-independent, but we will consider them being "weakly 
frequency dependent" in some sense. We then proceed by approximating their 
frequency dependence with "effective column densities". 



where Q is either Ei or D2, and Qqt is the respective quantity in the opticahy 
thin regime. Index a runs over the hst of species, which normally includes H I, 
He I, and Hell, but may also include other species, for example such as H2 
or dust. ^ Any rate of a specific physical process (such as photoionization 
rates, photoheating rates, etc) will now depend on those N^g^ . For the case of 
just three species (say, Q;=Hl,Hel,Hell), we precompute these rates with fine 

^ In practice we do a slightly more complicated trick, and assume that NeS = 
Ni + XN2 where X is the abundance of a species a (which can change during the 
hydrodynamic time-step), and are now fixed functions of position during the 
time-step (independent of time during the hydrodynamic time-step) . The quantity 
N2 is defined as N2 = nL, where n is the number density of a given species at a 
given point, and L is the size of the resolution element. The latter quantity is not 
defined precisely, but rather up to a factor of a few, but the numerical scheme is 
quite insensitive to its specific value: we varied it by a factor of 10 with little change 
to the final solution. 
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frequency dependence of Jy and g^^ store them in a 3D table, and then simply 
pick up a number from a table when it is needed. For four or more species 
rates need to be recalculated "on the fly" because the four-dimensional table 
is not practical. Tests show that precomputing gives about a factor of 10-30 
speed up. 

In order to use ansatz (14), we need to define the two quantities Ei^,^ and D2,v 
First, we note that since 

{Jv)q — JuEi,, 



then 

El, = Ei^y + y-D2,. 



The simplest way to split the "full" function E into Ei and D2 is by the 
following two equations (we put tildes over Ei and D2 because, albeit the 
simplest, this scheme does not work in practice): 



^adE^^, _ d ( 1 dE^,,,}^l 



c 




c dt dx^ \kj, dx^ 

adb2,u d ( 1 d]D2,vh'l^ in, - 



where /i*/^ = S^^ /3 and h2i, is the original Eddington tensor {h]^ from eq. [10]). 
Note, that the second equation now contains and not ki, because the Ji, 
dependence was factored out. 

However, the scheme (16) does not work, because £'2,^ = 9^02/ Jv can then be 

negative and positive, and numerical truncation errors will prevent the exact 
cancellation for E^ = -Ei^jy + £'2,1^ outside the ionization front, where E^, should 
be exactly zero. Therefore, we use a different splitting: 



c dt dx^ \k,y dx^ J ^ ^1 
adD2,v d ( 1 dD2,uh'lnu\ ^ ^ _ . 

Quantity is not specified here, it is only required to go to when the mean 
free path is much shorter than the size of the periodic computational box, and 
to 1 when the means free path is much larger than the box size. Irrespective of 
the choice of A,^, periodic boundary conditions of the computational box will 
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guarantee that our solution will deviate from the solution in the exact case 
of non-periodic universe when the mean free path is close to the size of the 
computational box. Therefore, the specific choice for is not that important. 
In this paper we adopt K^, = exp(— (kj^)L), where 

^ If cPxi 
Att J {x — xi)"^ 

(the latter integral is apparently independent of x). We verified that other 
choices of give different solutions when the mean free path is close to the 
box size, but converge to the same solution when the mean free path is much 
shorter or much larger than the box size. 

We emphasize that this procedure is only required in order to reduce the 
number of frequency bins in the calculation. It can be entirely circumvented 
if the full frequency space coverage can be afforded. 

The scheme presented above can be easily generalized on a case when there 
are several kinds of sources, each with its own spectral shape g^. 



3.3 Numerical implementation 



Equations (17) are then solved by the following simple semi-explicit scheme. 
Let Q be either Ei or D2 (we omit the frequency dependence in this section 
for the sake of clarity) , 



d (idQh 



dx^ \K dx^ J 

with the appropriate tensor h'^\ P is the source term {ipA for Ei or — p*A,^ 
for D2), and = cdt/{2a). Then equations (17) are shorthanded as 

^ = A2-«:g + p. (18) 

OK 

Let the superscript (n) label the n-th time step in the time-discretized version 
of equation (18). We then adopt a semi-explicit scheme for solving equation 
(18), in which we treat the Laplacian and the source term explicitly, and 
absorption implicitly (since the time scale for absorption can be extremely 
short) , 

/S(n+1) 
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which can be reduced to the exphcit expression for Q^'^'^^\ 

g("+i) = qW + 5qW, (19) 



and 

An) 



The time step 5^ is determined from the standard Courant condition for the 
diffusion equation, 

3 

5^ = Cmin(-Ax\), 



where Ax is the mesh cell size, and C — 0.9 is the Courant number. 

Our scheme is implemented in the Soften Lagrangian Hydrodynamics (SLH) 
code (Gnedin, 1995; Gnedin & Bertschinger, 1996; Gnedin, 2000). The SLH 
code follows all physical quantities on a moving deformed mesh specified by 
a coordinate transformation from the quasi-Lagrangian space (where the 
mesh is uniform) into the real space = x^{t,q^). In the quasi-Lagrangian 
space the Laplacian can be represented by using the deformation tensor 

A' = — 



its inverse Bf, and its determinant A: 

1 d 



Adq'' 



k + Ak ' ^ dq'- 



where we also introduced "artificial viscosity" An = e/Ax to prevent the 
denominator from becoming zero. The results are insensitive to the choice of 
e as long as it is less than about 0.03. The cell size Ax is now defined as 

Ax"" ^2>(B^,B\5ij5^^y\ 

The description above completes our numerical scheme except for the choice of 
the "reduced speed of light" c. We determine it from the following condition: 

c = c X min(l, N^^St/At) (20) 
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where At is the hydrodynamic time-step during which equations (17) are 
solved, 5t is the time-step for equations (17) determined from the Courant 
condition, and N^^y^ is the upper hmit on the number of time-steps in one 
hydro time-step. This is equivalent to solving equations (17) in N-^^ time 
steps at most - if more time steps is required, the integration is interrupted 
at the moment when the maximiim allowed number of time steps is reached. 
Our scheme is fairly insensitive to the choice of A^max, as discussed in the next 
subsection. 

Thus, the whole scheme for solving RT within the SLH code is the following. 
At each hydrodynamic time step: 

• gravity and hydro are updated. 

• equations (17) are solved for A^max steps at four frequency bins (i/=0,H I, He I, He ll) . 
Then n'^^ are calculated. 

• equations for gas temperature and ionization state are solved given n'^ . 



3.4 Signal propagation speed 



As we discussed in §3.1, the RSL approximation by itself does not limit the 
signal propagation speed. However, an adopted limit on the number of time- 
steps A'lnax does limit the speed with which a signal can propagate across 
the mesh. Namely, since we are using a 3-point finite differencing scheme for 
computing the Laplacian, in one time-step a signal can only propagate one cell 
on our quasi-Lagrangian mesh. Let us first assume that our mesh is uniform in 
real space as well, i.e. that our quasi-Lagrangian space coincides with the real 
space. In that case during one hydrodynamic time-step At the signal can only 
travel the physical distance Ar = aLiVmax/A^box- where L is the comoving size 
of the computational box, and A'box is the number of cells along one dimension 
of the box. Thus, the propagation speed of a signal over the box is limited to 

or 

'^max O-^sA^max f21) 



where L3 = L/3Mpc and Atio = At/10 Myr. For example, for the test shown 
in Fig. 1 this number is about 1.2A^jnax- The test was performed with A^max = 1 
and A^max = 4 and we found no noticeable difference between the two cases, 
which is not surprising since the top speed of the ionization front in this case 
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is only about 4000 km/s. This speed is however much larger than A^max times 
the sound speed in the gas (16km/s behind the ionization front and 4km/s 
outside it), which illustrates that the RSL approximation does not limit the 
propagation speed of the signal. 

Thus, one can in principle adjust iVmax in equation (21) so that the signal 
propagation speed is always fixed to c (unless it would require A^max < !)• 
However, for realistic cosmological simulations presented in §4.5 we do not do 
that, because the SLH code allows for severe deformations of the computa- 
tional mesh in real space, and in this case the signal propagation speed is only 
fixed in the quasi-Lagrangian space for a fixed value of A^maxi but varies at 
different spatial locations in real space. Instead, we keep A^max fixed through- 
out the simulation, and we have tried two values, iVmax = 20 and iVmax = 100, 
which give visually indistinguishable results. Incidentally, during the whole 
simulation with A^max = 20 the signal propagation speed Vmax never deviated 
from c by more than a factor of 10 (and it is respectively 5 times larger for 

A^max - 100). 



3. 5 The choice of the Eddington tensor 



The scheme presented in the previous sections is complete given the Eddington 
tensor k^J. The devil of course is in the choice for the Eddington tensor, because 
it allows to close the moment hierarchy. 

The Eddington tensor can be written in a closed form as 

pij 



where 

P'J = I d'x,Mx,)e-^-^^'''^^ '''^ , (23) 



where t^{x,Xi) is the optical depth between points x and Xi. Evaluation of 
this integral requires 0{N^^^) operations, and thus this form is unsatisfactory. 
Thus, we need to make a further approximation to achieve a favorable 0{N) 
scaling. However, this approximation has to be introduced carefully so that 
the physical meaning of equation (23) is not violated. Specifically, equation 
(23) describes the fact that radiation propagates from a set of discrete sources 
along straight lines, and any new approximation has to preserve this property. 

In this paper we propose to use Eddington tensors calculated in the optically 
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thin regime, 



P« = /rf'..P.(fO^^^^|%^ (24) 

(in this case P*-' and become frequency independent). The latter integral 
can now be computed in 0{N) operations by standard algorithms used to 
compute high-resolution gravitational forces such as P^M, Tree, or Adaptive 
Mesh Refinement. 

It is important to underscore that in our approximation the physical meaning 
of equation (23) is preserved - radiation still propagates from a set of discrete 
sources along straight lines - and only the relative weight of different sources is 
misestimated. The latter however does not mean that the number of ionizing 
photons emitted by a given source is misestimated as well. By separating 
the computation of the Eddington tensor from the first two moments of the 
transfer equation, we in effect separate two physical processes: absorption and 
emission of photons (which is described by the two moments) and propagation 
of photons along the straight lines from the source locations (which is described 
by the Eddington tensor). The first we compute almost exactly (subject to the 
RSL approximation and finite number of frequency bins) , while the second we 
compute in the optically thin regime. In other words, we insure the exact 
conservation of the number and the flux of photons, but may occasionally 
advect some flux in the wrong direction. 

We also stress that out approximation differs significantly from other schemes, 
like the diffusion approximation. In the latter the Eddington tensor is com- 
puted locally, and it thus misses the non-locality present in equations (23) 
and (24). For example, the diffusion approximation will grossly fail for the 
case when a small Hll region is overrun by a larger Hll region from a more 
distant but a much stronger source, whereas our approach will give a quite 
accurate representation of the radiation field in that case (this test is discussed 
in the next section). 



4 Testing the scheme 

4-1 Single source tests 

In order to verify the accuracy of our method, we have compared the results 
obtained with the 3D code to the predictions of the exact solution in spherical 
symmetry. 
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In spherical symmetry equation (7) with no time derivative term and for a 
point source at the center {ipi, — has an exact solution, 



fuit,r,n) 



r 

J k,{t,r')dr' 



S{n — Cr 



Our implementation of a spherically symmetric code has 1000 points per 
decade in the radial direction (which has been verified to be sufficient to 
converge to 1% accuracy) and exactly the same frequency dependence and 
ionization and heating rates as the 3D SLH code. 

Figure 1 now illustrates how our scheme works for a case of a point source in a 
uniform medium. In this test the point source is located at the center of a 64^ 
mesh in a 70 kpc box with periodic boundary conditions. The source spectrum 
contains only hydrogen ionizing photons. We notice that the approximate 
3D solution converges to the exact solution in about 10 hydrodynamic time 
steps, after which time the 3D scheme reproduces accurately not only the 
position of the ionization front, but also its detailed temperature and ionized 
fraction profiles. The deviation at the initial moment is most likely due to 
finite resolution of the 3D code. At later times the ionization front leaves the 
computational box and the solution properly approaches the optically thin 
regime (open symbols in Fig. 1) when most of the box has a uniform radiation 
background and a proximity region exists around the source. 

Figure 2 shows in two panels the same test with 100 and 10,000 weaker source 
respectively. Again, not only the position, but also the detailed profile of the 
ionization front is reproduced. In the latter case the source is so weak (by 
construction) that it reaches its Stromgen sphere before the ionization front 
leaves the box. We continued the weak source test for another factor of three 
in time and the solution remains undistinguishable from the last line shown 
with no trace of unphysical numerical diffusion propagating outside of the 
Stromgen sphere. 

Figure 3a shows the effect of changing the hydrodynamic time step. Three sets 
of filled symbols correspond to the hydrodynamic time step changed by a factor 
of 1/2, 1, and 2 relative to Fig. 1. We again notice that about 10 hydrodynamic 
time steps after switching on of the source are sufficient to achieve complete 
convergence. Figure 3b also tests the performance of our scheme in the case 
when the ionizing spectrum extends to higher energies, thus being able to 
ionize He I and Hell as well as Hi. In that case the agreement for helium 
is somewhat worse than for hydrogen, and it takes about 30 hydrodynamic 
time steps to converge on the helium profiles. This is however is still very 
small compared to the total number of hydrodynamic time steps in a realistic 
simulation, so that the compound error is expected to be negligible. 
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Fig. 1. Propagation of the ionization front into the uniform neutral medium with 
n = 2.5 X 10~^ cm~^ from a point source with iVph = 10^^ s~^. Sohd hues show the 
exact spherically symmetric solution, whereas grey symbols show the 3D approxi- 
mate solution after 3, 6, 12, 25, 50, 100, and 200 hydrodynamic time steps (in the 
order of the increasing darkness of a symbol). Two sets of open symbols show the 
3D solution after 700 and 1000 hydrodynamic time steps. The latter moment corre- 
sponds to the time when the entire computational box is optically thin. The upper 
panel shows the gas temperature, the middle panels shows the hydrogen ioniza- 
tion fraction on a linear scale, and the bottom panel shows the hydrogen ionization 
fraction on a logarithmic scale. 



Figure 4 illustrates another two tests in which the density distribution is not 
uniform but falls off away from the source as the first or the second power of 
radius. In the latter case the profile of the ionization front at earlier times is 
not well reproduced due to the finite resolution of the 3D scheme. 
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Fig. 2. The same as Fig. 1 except for iVph = 10^° s"^ (a) and iVph = 10"^^ s"^ (6). In 
this case only first 200 hydrodynamic time steps are shown. 
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Fig. 3. The same as Fig. 1 except ( a) three sets of solutions shown corresponding to 
three choices for the hydrodynamic time step; (b) the source spectrum is a power-law 
with the slope of -1.25 and thus includes Hei and Hell ionizing photons. In the 
latter case the middle panel shows the He I ionization fraction (filled symbols). Open 
symbols mark He II ionization fraction at the last output. Only 6th, 25th, and 100th 
hydrodynamic time steps are shown. 



4-2 Shadowing 



The tests from the previous section demonstrate that our scheme performs rea- 
sonably well for a single source in a spherically symmetric case. However, the 
choice of the optically thin Eddington tensor is, obviously, an approximation, 
and wc can identify regimes where it is likely to fail. This is best illustrated 
by Figure 5, which shows the two eigenvectors of the optically thin Eddington 
tensor for a single source. The box has a unit side, and is located so that its 
center coincides with the point (0,0,0). The source is located at the center. We 
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Fig. 4. The same as Fig. 1 except that the density is not uniform by falls off as r ^ 
(a) or (6) with a small core at the center. 

show the X — y plane of the box with two eigenvectors laying in that plane. 
The eigenvectors are scaled by the square root of their respective eigenvalues 
(due to symmetry, the third eigenvector is directed along z direction, but the 
respective eigenvalue is zero). Their signs therefore are undetermined, i.e. each 
vector can be turned around 180 degrees. 



In this test case it is easy to find the exact Eddington tensor - it has only 

one nonvanishing eigenvalue (therefore equal unity everywhere because the 
trace of the Eddington tensor is identically equal unity) and the corresponding 
eigenvector at every point in space points towards (or away from) the source 
(i.e. h^^ = e].ei). The grey vectors in Fig. 5 show this eigenvector. However, in 
the optically thin case there exist another eigenvector, shown in black in Fig. 5. 
This eigenvector arises due to periodic boundary conditions as a contribution 
from periodic images of the source. As the result, the second eigenvector will 
allow for the propagation of the ionization front in the direction tangential to 
the source (albeit with a lower speed). 

Figure 6 serves to illustrate this point. Its lower left panel shows the one 
quadrant in the midplane from a test shown in Fig. 1. Dark color shows the 
ionized region, and light color shows the still neutral medium. The lower right 
panel shows the same test (the panel is flipped around the y axis) but with 
the 1000 times denser wedge inserted on the way of the ionization front. The 
wedge is supposed to shield the medium behind it from the ionizing radiation, 
however we can see the ionization front propagating behind the wedge with 
about 1/5 of the normal ionization front speed. This propagation is caused by 
the second eigenvector of the Eddington tensor, and is therefore unphysical 
- it is a manifestation of the optically thin Eddington tensor calculated in a 
periodic universe. 

To verify that our conclusion is indeed correct, we have performed the same 
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Fig. 5. Two eigenvectors of the optically thin Eddington tensor in one quadrant of 
2 = section of the computational box. The source location is labeled by a filled 
circle. 

test in a twice larger box (128^ instead of 64^) shown on the upper left panel 
of Fig. 6. In that case all the periodic images are twice further away, the 
contribution from them to the Eddington tensor is 4 times smaller, and we 
expect the tangential ionization front to propagate with (1/4)^/'^ = 0.6 times 
lower speed, as indeed can be observed from the figure. 

Finally, we show the same test but calculated with the exact Eddington tensor 
on the upper left panel. In this case the ionization front does not propagate 
behind the wedge except by a very weak numerical diffusion (which is not 
important in a realistic simulation). 

These tests illustrate the main shortcoming of the proposed scheme: the slow 
propagation of the ionization fronts behind the shadow. We emphasize that 
this propagation is quite slow, significantly slower than the ionization front 
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Fig. 6. The distribution of neutral hydrogen in the one quadrant of the midplane of 
the computational box for four different tests: uniform medium {lower left panel), 
uniform medium with a 1000 times denser wedge (lower right panel), uniform 
medium with a wedge in a twice larger (128^) box (upper right panel), uniform 
medium with a wedge with the exact Eddington tensor (upper left panel). All pan- 
els arc flipped so that the source is located at the center of the plot. The side 
bars give the correspondence between the color and the decimal log of the neutral 
hydrogen density. 

speed. We notice this effect only because there should be no propagation of 
the radiation behind the shadow whatsoever. Our tests are not sufficient to 
evaluate whether this is a significant drawback or it is not important in realistic 
simulations - we expect that future work will shed light on this question. 

In the proposed framework, we are not aware of the acceptable way to circum- 
vent this deficiency of the scheme. We have tried to experiment with different 
choices of the Green functions for calculating the Eddington tensor, but found 
no noticeable improvement - with the following simple reason: if we adopt a 
Green function which falls faster than (a short-range one) in order to 
reduce the contributions from the distant periodic images, we at same time 
introduce deformations into the Eddington tensor inside the computational 
box, because of the periodicity of boundary conditions. These deformations 
may lead to even larger errors than the original optically thin Eddington ten- 
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sor. Significant improvement over the proposed scheme is only possible if the 
Eddington tensor is estimated more accurately, for example, by calculating 
the exact tensor at a number of selected points inside the computational box. 
However, it seems unlikely this could be done with acceptable computational 
cost. 

We again emphasize here that the fact that the Eddington tensor is calculated 
in the optically thin regime does not compromise shadowing by itself, but only 
due to imposed periodic boundary conditions. The reason that the shadowing 
works with the optically thin Eddington tensor is because the absorption of 
ionizing photons by the obstacle is computed essentially exactly (up to the 
RSL approximation, which becomes exact in the limit of a very dense obstacle). 
The optically thin approximation does make an error by assigning the incorrect 
values for the Eddington tensor behind the obstacle, but because the moments 
equations insure that the ionizing photon number density behind the obstacle 
is zero, the error in the Eddington tensor becomes irrelevant. 

4.3 Multiple sources 

It may appear that the choice of the optically thin Eddington tensor will also 
lead to a serious error in a case of several sources inside the computational box. 
Indeed, let's consider two sources inside a computational box: a strong one and 
a weak one, and let's assume that the two sources switch on simultaneously. 
Until their respective ionization fronts meet, each will have a spherical Hll 
region around it, and inside each of the Hll regions the Eddington tensor will 
be of the form n^n^ where n is the unit vector pointing toward the respective 
source. However, the proposed scheme adopts an optically thin Eddington 
tensor, which means that at each point in space the Eddington tensor is given 
by equation (22) with 

where fi and f2 are coordinates of the two sources in a reference frame tied 
to a given point in space, and Si and S2 are their respective luminosities. We 
assume that Si > 82- 

Near the strong source the contribution from the weak source is small, and 
thus the optically thin Eddington tensor is a good approximation to the exact 
one. However, near the weak (second) source the contribution from the first 
source could be important: even if r2 < ri. Si > S2. However, the physics 
of the ionization fronts is such that this is not a significant problem. We can 
illustrate that with the following simple estimate. 
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Fig. 7. The distribution of neutral hydrogen in one quadrant of the midplane of 
the computational box for four different cases. The upper row shows the two source 
test with the weak source located at (-0.5,-0.5,0) and in the bottom row the weak 
sources is at (-0.25,-0.25,0). The strong source is at (0,0,0), and the ratio of source 
luminosities is 10. Left columns show the case with the exact Eddington tensors 
(which can only be calculated easily before the two ionization fronts overlap), and 
right columns show the test with optically thin Eddington tensors. In both cases 
the difference is noticeable but not large. 

The radius of the ionization front scales as a cubic root of the source luminosity. 
Thus, at the moment when the ionization fronts from the two sources overlap, 
the ratio of the front radii is 

RI/R2 = {Sl/S2y/\ 

The contribution to the tensors from each of the sources is proportional to 
S/r'^, i.e. the ratio of the contributions from the weak and the strong sources 
to the Eddington tensor within the weak source's H ll region is 

hi' _ S^rl 
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Fig. 8. The distribution of neutral hydrogen in one quadrant of the midplane of the 
computational box for the two-source test with optically thin Eddington tensors at 
four different times. 

This contribution is greater than 1 (i.e. the weak source dominates its own 
Eddington tensor in the optically thin regime regime) when 

\ K2 J \ tX\ J Ki 

or 

i?2 Ri \SiJ 

Even at the time of overlap, when ri ~ Ri, most of the volume inside the 
HII region around the weak source is dominated by the weak source, because 
1/6 is too weak a power law. Before the overlap, ri > Ri (ri is the distance 
from the strong source to a point inside the weak source H ll region) , so the 
constraint is even weaker. 

Figure 7 illustrates this conclusion. The two two-source tests presented in the 
figure show that the effect of the stronger source on the propagation of the 
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Fig. 9. The number of ionizing photons per baryon as a function of time for the 
optically thin case {dashed line) and for the full calculation [solid line) for the test 
presented in Fig. 1. The dotted line shows the optically thin case minus one (one 
photon per baryon is used to ionize the box). 

ionization front around a weaker source are small albeit still noticeable. We 
also show in Figure 8 consequent evolution for this test for the approximate 
solution. We are unable to compute the exact solution in this case because we 
cannot calculate the exact Eddington tensor when two H ll regions overlap. 



4-4 Optically thin regime 



One of the attractive features of this scheme is that the transition to the 
optically thin regime (or, more precisely, to the regime when the mean free 
path is much larger than the size of the computational box) is automatic, 
whereas in many alternative approaches it is not so trivial. Fig. 1 shows this 
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regime for one of our tests with open symbols. In Figure 9 we show the ionizing 
background (i.e. the volume averaged specific intensity) measured in units of 
ionizing photons per baryon in the computational box for the optically thin 
case (no absorptions) and for one of our single source tests. The dotted hne 
in this figure shows the optically thin case minus one, which approaches the 
full solution at later times. This illustrates the conservation of photons as 
the solution switches from the optically thick to the optically thin regime (in 
this test, recombinations are not important, and only one ionizing photon per 
baryon is needed to ionize the whole box). 

This test illustrates that our scheme automatically insures the exact con- 
servation of the number of photons as the mean free path crosses the box 
boundaries. 



4.5 Realistic simulations 

While the main purpose of this paper is to present the new and accurate ap- 
proximation for modeling continuum radiative transfer, we give here also an 
example of how it performs for the problem it was developed for, cosmolog- 
ical reionization. For this purpose we have repeated one of the simulations 
from Gnedin (2000), namely the run labeled "N64X2_A". That simulation 
contained 64^ dark matter particles, an equal number of baryonic cells and, 
towards the end of the calculation, roughly the same number of stellar par- 
ticles in a computational box with the size of 2h~^ comoving Mpc. A reader 
can find a full description of this simulation in Gnedin (2000) and we do not 
repeat it here for brevity. This test allows us not only to check the performance 
of our Optically Think Variable Eddington Tensor (OTVET) approximation 
in a realistic cosmological simulation, but also to compare it with the much 
less accurate Local Optical Depth (LOD) approximation. The OTVET ra- 
diative transfer simulation presented here consumed roughly three times the 
computational time than the one using LOD transfer. 

In Figure 10 we show the evolution of the volume averaged ionizing intensity 
for two simulations. The two simulations have identical parameters and the 
same physical content - the only difference is in the method used to solve for the 
radiative transfer. We notice here that the two simulations are qualitatively 
similar, and reionization (the sharp rise in the specific intensity at around 
z — 7) occurs almost at the same time in both simulations. This should be 
expected because the LOD approximation is designed to give the correct time 
for the overlap of H II regions, which marks the time of reionization. However, 
there are also noticeable differences, which we attribute to the inaccuracy of 
the LOD approximation. 
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Fig. 10. The evolution of the spatially averaged ionizing intensity Jjj j (measured in 
conventional units of 10"^''^ erg/ cm^/ sec/ Hz/ rad) in two simulations of a represen- 
tative cosmological model: using the Local Optical Depth (LOD) approximation of 
Gnedin (2000, dotted line ) and Optically Thin Variable Eddington Tensor (OTVET; 
solid line) approximation developed in this paper. 

In order to investigate these differences further, we show in Figure 11 shoes 
from two simulations at three different times - before the overlap of H II regions, 
during the overlap, and after the overlap. We notice that while large-scale fea- 
tures are similar in two simulations, small-scale structure is quite different. 
The main difference is the excess of low density neutral gas in the LOD ap- 
proximation. This is indeed expected, because the LOD approximation makes 
the largest error in the regions with the optical depth is about unity. 

It appears from Fig. 11 that this error is systematic - neutral regions are always 
larger in the LOD approximation - and thus the LOD approximation performs 
worse than one of the authors (NG) hoped for. On the other hand, the LOD 
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Fig. 11. The distribution of neutral hydrogen in a reahstic simulation of cosmological 
rcionization at three different times as shown in three rows: before the overlap of H II 
regions (top row), during the overlap (middle row), and after the overlap (bottom 
row). The left column shows the highly approximate solution with the Local Optical 
Depth (LOD) approximation from Gnedin (2000), and the right column shows the 
present scheme (OTVET). 

approximation only reproduces the ionization fronts speed on average, i.e. it 
does not guarantee a correct speed for every ionization front. The fact that 
the large-scale structure of ionization fronts agrees quite well between the two 
simulations indicate that the ionization front speed in the LOD approximation 
is calculated more accurately than could be expected a priori. This observation 
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would allow to improve the LOD approximation by correcting the systematic 
deviation on the r ~ 1 regime. 

Another significant difference between the two approximations is a higher 
asymmetry of Hll regions at early times in the LOD approximation. This 
difference arises because the LOD approximation estimates the total opti- 
cal depth to the source from its local density. Somewhat denser material 
becomes therefore artificially more difficult to ionize. Consequently, I-fronts 
wrap around dense material unphysically strongly giving rise to the strongly 
distorted ionization fronts observed in Figure 11. 

Finally, at later times (after the overlap) the two solutions are quite similar 
(except for small errors in the LOD approximation in the r ~ 1 regime) 
because both approximations become exact in the optically thin regime. 



5 Conclusions 

We present a new approach to numerically model continuum radiative transfer 
based on the Optically Thin Variable Eddington Tensor (OTVET) approxi- 
mation. Our method insures the exact conservation of the photon number and 
flux (in the explicit formulation) and automatically switches from the optically 
thick to the optically thin regime. It scales as N log N with the number of hy- 
drodynamic resolution elements and is independent of the number of sources 
of ionizing radiation (i.e. works equally fast for an arbitrary source function). 

Wc also describe an implementation of the algorithm for the Soften Lagrangian 
Hydrodynamic code (SLH). We present extensive tests of our method for sin- 
gle and multiple sources in homogeneous and inhomogeneous density distri- 
butions, as well as a realistic simulation of cosmological reionization. 

The method fully solves the cosmological radiative transfer equation in a con- 
servative fashion via moment equations. The only approximation is our use of 
optically thin Eddington tensors. One may envision to employ a Monte Carlo 

or short characteristic methods to supply more accurate tensors for applica- 
tions in which an even more accurate knowledge of the isotropy of radiation 
fleld is required. 

The main shortcoming introduced by optically thin Eddington tensors with 
periodic boundary conditions is the unphysical propagation of ionization fronts 
behind an optically thick obstacle. This diffusive effect is however slow and 
negligible for the typical R-type ionization fronts of interest in numerical 
cosmology. 
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A realistic simulation of cosmological reionization with our current implemen- 
tation of the new method took approximately three times the computing time 
than a simulation without radiative transfer. However, there are several op- 
timizations which we have not yet implemented, which will shorten this time 
by at least a factor of two. 
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